Introduction
Packages
library(tidyverse)
library(plotly)
library(scales)Input Parameters
S <- 1370
A <- 204
B <- 2.17
K <- 3.86
ai <- 0.62
ab <- 0.25
Temperature <- rep(0,300)
Temp <- rep(0,100)
Temp2 <- rep(0,180)
Temp5 <- rep(0,300)
J <- rep(0,100)
Albedo <- rep(0,300)
t <- c(1:300)
Zones <- seq(-89, 89, by = 2)Fit from Textbook Spreadsheet
zones <- c(85,75,65,55,45,35,25,15,6)
coszones <- cospi(zones/180)
sunWt <- c(0.5,0.531,0.624,0.77,0.892,1.021,1.12,1.189,1.219)
df <- data.frame(zones,sunWt)
f <- function(zones,d,n,k){
d*cos(n*zones)^2+k } #cos or cos^2 ??
Fit <- nls(sunWt ~ f(zones,d,n,k),data= df, start=list(d=1,n=0.015,k=0.3))
ggplot(df, aes(zones,sunWt))+geom_point()+geom_smooth(aes(y=f(zones,0.7768699,0.0164348,0.4617747)))+xlab("Latitude")+ylab("Solar Weight Factor")
Functions
Func <- function(x){
0.7768699*cos(0.0164348*x)^2+0.4617747
}
gauss <- function(x,m,sd,b){
((24+b)/(0.00798*sqrt(2*pi*sd^2)))*exp(-(x-m)^2/(2*sd^2))-b
}
Sun1 <- function(x,a){a*x/100}
Sun2 <- function(x){1370*(sinpi((x+90)/180))^2}
Sun3 <- function(x){1370-(((1370)/(sqrt(2*pi*1^2)))*exp(-(x-50)^2/(2*1^2)))}
Sun4 <- function(x){ifelse(x==50,1370/3,1370)}
Sun5 <- function(x){(1/100) * (abs(x-150)+25)}
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
Incident <- function(x,y){x*y/4}
Step <- function(x,c){ifelse(x<c, 0.6, 0.3)}
alb <- function(x,a,b){
(-exp(2.2*(x+10/2.2))/(exp(2.2*(x+10/2.2))+1))*(a-b)+a}\[ S_1(t)= \frac{1370}{100}t \] \[ S_2(t)= 1370(\sin^2(t+90)) \] \[ S_3(t)= 1370\Big(1-\frac{1}{\sqrt{2\pi}}e^{-(t-50)^2/2}\Big) \] \[ S_4(t)=1370\Big(1-\frac{1}{3}\delta(t-50)\Big) \] \[ S_5(t)=\frac{1}{100}(|t-150|+25) \] \[ S_6(t)=1370\Big(1+\frac{1}{10} cos(t)\Big) \] \[ a(T)=\frac{e^{2.2T+10}}{e^{2.2T+10}+1}(b-a)+a \]
ggplot(data.frame(x= seq(-15,15, by=0.1),y=alb(seq(-15,15, by=0.1),ai,ab)))+geom_line(aes(x,y),colour="blue")+xlab("Temperature")+ylab("Albedo")
Ti <- gauss(Zones,0,50,31.6)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
Rin <- Incident(S,SunWt)
T <- gauss(Zones,0,50,31.6)
a <- alb(T,ai,ab)Run 0
for(i in c(1:300)) {Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T,ai,ab)
Temperature[i] <- Tm
}
data1 <- data.frame(Zones,T,a,Ti)
plot1 <- ggplot(data1,aes(Zones)) +
geom_line(aes(y=T, colour = "Temperature",linetype="Temperature")) +
geom_line(aes(y=a*25, colour = "Albedo",linetype="Albedo"))+
geom_line(aes(y=Ti, colour = "Initial Temperature",linetype="Initial Temperature"))+xlab("Latitude")+ylab("Temperature")+scale_colour_manual(name="Legend",values=c("blue","red","red"))+scale_linetype_manual(name="Legend",values=c("Temperature"=1, "Albedo"=1, "Initial Temperature"=2))
plot1
data2 <- data.frame(t,Temperature)
plot2 <- ggplot(data2) +
geom_line(aes(t, Temperature),colour = 'green')+xlab("Time")+ylab("Mean Temperature")
plot2
Run 1
TEMP1 <- matrix(NA, nrow=90, ncol=100)
for(j in c(1:100)){
T <- gauss(Zones,0,50,31.6)
a <- alb(T,ai,ab)
S <- Sun1(1370,j)
Rin <- Incident(S,SunWt)
for(i in c(1:300))
{Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T,ai,ab)}
TEMP1[,j] <- T
J[j] <- j
Temp[j] <- Tm
}
data3 <- data.frame(J,Temp)
plot3 <- ggplot(data3,aes(J,Temp))+geom_line(aes(J, Temp),colour = 'red')+xlab("Time")+ylab("Mean Temperature")+ggtitle("Non Dynamic Mean Temperature (T(t) independent from T(t-1))")
plot3
plot_ly(z=~TEMP1)%>% add_surface() %>% layout(
title = "With Initialization", scene = list(
xaxis = list(title = "S/100"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "T") )) Run 2
T2 <- gauss(Zones,0,50,31.6)
TEMP2 <- matrix(NA, nrow=90, ncol=180)
a <- alb(T2,ai,ab)
Sarr <- rep(0,180)
for(j in c(1:180)){
S <- Sun2(j)
Rin <- Incident(S,SunWt)
for(i in c(1:300))
{Tcos <- cosZones*T2
Tm <- sum(Tcos)/sum(cosZones)
T2 <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T2,ai,ab)
}
Sarr[j] <- Sun2(j)
TEMP2[,j] <- T2
J[j] <- j
Temp2[j] <- Tm
}
data3 <- data.frame(J,Temp2)
plot3 <- ggplot(data3,aes(J,Temp2))+geom_line(aes(J, Temp2),colour = 'red')+xlab("Time")+ylab("Mean Temperature")+ggtitle("Sinusoidal Perturbation")
plot3
plot_ly(z=~TEMP2) %>% add_surface() %>% layout(
title = "Without Initialization",scene = list(
xaxis = list(title = "S/100"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "T") ))Run 3 & 4
TEMP3 <- matrix(NA, nrow=90, ncol=200)
T3 <- gauss(Zones,0,50,31.6)
a <- alb(T3,ai,ab)
for(j in c(1:200)){
S <- Sun3(j)
Rin <- Incident(S,SunWt)
for(i in c(1:300))
{Tcos <- cosZones*T3
Tm <- sum(Tcos)/sum(cosZones)
T3 <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T3,ai,ab)
}
TEMP3[,j] <- T3
J[j] <- j
Temp[j] <- Tm
}
data3 <- data.frame(J,Temp)
plot3 <- ggplot(data3,aes(J,Temp))+geom_line(aes(J, Temp),colour = 'red')+xlab("Time")+ylab("Mean Temperature")+ggtitle("Gaussian Perturbation")
plot3
#Run4: senza inizializzazione di (T,a)
TEMP4 <- matrix(NA, nrow=90, ncol=200)
T4 <- gauss(Zones,0,50,31.6)
a <- alb(T4,ai,ab)
for(j in c(1:200)){
S <- Sun4(j)
Rin <- Incident(S,SunWt)
for(i in c(1:300))
{Tcos <- cosZones*T4
Tm <- sum(Tcos)/sum(cosZones)
T4 <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T4,ai,ab)
}
TEMP4[,j] <- T4
J[j] <- j
Temp[j] <- Tm
}
data3 <- data.frame(J,Temp)
plot3 <- ggplot(data3,aes(J,Temp))+geom_line(aes(J, Temp),colour = 'red')+xlab("Time")+ylab("Mean Temperature")+ggtitle("Delta Perturbation")
plot3
plot_ly(z=~TEMP3) %>% add_surface()%>% layout(
title = "Without Initialization", scene = list(
xaxis = list(title = "S/100"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "T") ))plot_ly(z=~TEMP4) %>% add_surface()%>% layout(
title = "Without Initialization",scene = list(
xaxis = list(title = "S/100"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "T") ))Hysteresis Cycles
data <- data.frame(Sarr,Temp2)
plot <- ggplot(data,aes(Sarr,Temp2))+geom_point(aes(Sarr, Temp2),colour = 'yellow')+ylab("Mean Temperature")+xlab("S_2(t)")+ggtitle("Mean temperature vs. Sinusoidally Fluctuating Incoming Radiation")
plot
Temp5 <- rep(0,300)
T5 <- gauss(Zones,0,50,31.6)
a <- alb(T5,ai,ab)
Sarr <- rep(0,300)
for(j in c(0:300)){
S <- Sun5(j)
Rin <- 1370*Incident(S,SunWt)
for(i in c(1:60))
{Tcos <- cosZones*T5
Tm <- sum(Tcos)/sum(cosZones)
T5 <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T5,ai,ab)
}
Sarr[j] <- Sun5(j)
J[j] <- j
Temp5[j] <- Tm
}
dataw <- data.frame(Sarr,Temp5)
plot5 <- ggplot(dataw,aes(Sarr,Temp5,group=1))+geom_point(aes(Sarr, Temp5),colour = 'yellow')+geom_segment(aes(x = Sarr[89], y = Temp5[89], xend = Sarr[90], yend = Temp5[90]),colour = 'yellow')+geom_segment(aes(x = Sarr[258], y = Temp5[258], xend = Sarr[259], yend = Temp5[259]),colour = 'yellow')+ylab("Mean Temperature")+xlab("S_5(t)")+ggtitle("Mean temperature vs. Linearly Fluctuating Incoming Radiation")
plot5
Embedding Daisies in EBM
#Input Parameters
S <- 1370
A <- 204
B <- 2.17
K <- 3.86
ai <- 0.62
ab <- 0.25
aW <- 0.75
aB <- 0.25
c <- 7
k <- 0.003265*0.75
T0 <- 20
D <- 0.3 # Death Rate
Zones <- seq(-89, 89, by = 2)
#Functions
Func <- function(x){
0.7768699*cos(0.0164348*x)^2+0.4617747
}
gauss <- function(x,m,sd,b){
((24+b)/(0.00798*sqrt(2*pi*sd^2)))*exp(-(x-m)^2/(2*sd^2))-b
}
Incident <- function(x,y){x*y/4}
alb <- function(x,a,b){ (-exp(1.6*x+10)/(exp(1.6*x+10)+1))*(a-b)+a} #naked
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
#Initialization
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
Rin <- Incident(S,SunWt)
T <- gauss(Zones,0,50,31.6)-6
w <- rep(0.5,length(Zones))
b <- rep(0.2,length(Zones))
u <- rep(0.3,length(Zones))
a <- w*aW+b*aB+u*alb(T,ai,ab)
Barr <- rep(0,500)
Warr <- rep(0,500)
Uarr <- rep(0,500)
Tarr <- rep(0,500)
I <- rep(0,500)
TEMP <- matrix(NA, nrow=90, ncol=500)
for(i in c(1:500)) {
S <- Sun6(i) # oppure costante S <- 1370
Rin <- Incident(S,SunWt)
Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
TEMP[,i] <- T
a <- alb(T,0.62,0.25)
I[i] <- i
Tarr[i] <- T[45]
}
plot <- ggplot(data.frame(Zones,w,b,u,T),aes(Zones))+geom_line(aes(y=w, colour = "% White"))+geom_line(aes(y=b, colour = "% Black"))+geom_line(aes(y=u, colour="% Bare Ground"))+geom_line(aes(y=T/5, colour="Temperature"))+ylab("Temperature & Albedo")+xlab("Latitude")+ggtitle("Without Daisies")+scale_colour_manual(name="Legend",values=c("brown","black","white", "red"))
plot
plot2 <- ggplot(data.frame(I,Barr,Warr,Uarr,Tarr),aes(I))+geom_line(aes(y=Barr,colour="% Black"))+ geom_line(aes(y=Warr,colour="% White"))+ geom_line(aes(y=Uarr, colour="% Bare Ground"))+geom_line(aes(y=Tarr/25,colour="Temperature"))+ylab("Temperature & Albedo")+xlab("Time")+ggtitle("Without Daisies")+scale_colour_manual(name="Legend",values=c("brown","black","white", "red"))
plot2
plot_ly(z=~TEMP)%>% add_surface() %>% layout( title="Without Daisies",
scene = list(
xaxis = list(title = "Time"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "T") )) a <- w*aW+b*aB+u*alb(T,0.62,0.25)
T <- gauss(Zones,0,50,31.6)-6
TEMP <- matrix(NA, nrow=90, ncol=500)
for(i in c(1:500)) {
S <- Sun6(i) # oppure costante S <- 1370
Rin <- Incident(S,SunWt)
Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
TEMP[,i] <- T
Tw <- T+c*(a-aW)
Tb <- T+c*(a-aB)
Fw <- 1-k*(T0-Tw)^2
Fb <- 1-k*(T0-Tb)^2
for(j in c(1:length(Zones))){
if(Fw[j]<0){Fw[j]=0}
if(Fb[j]<0){Fb[j]=0} }
w <- w+w*(u*Fw-D)
b <- b+b*(u*Fb-D)
for(j in c(1:length(Zones))){
if(w[j]<0.001){w[j]=0.001}
if(b[j]<0.001){b[j]=0.001} }
u <- 1-w-b
a <- w*aW+b*aB+u*alb(T,0.62,0.25)
Barr[i] <- b[45]
Warr[i] <- w[45]
Uarr[i] <- u[45]
I[i] <- i
Tarr[i] <- T[45]
}
plot <- ggplot(data.frame(Zones,w,b,u,T),aes(Zones))+geom_line(aes(y=w, colour = "% White"))+geom_line(aes(y=b, colour = "% Black"))+geom_line(aes(y=u, colour="% Bare Ground"))+geom_line(aes(y=T/5, colour="Temperature"))+ylab("Temperature & Albedo")+xlab("Latitude")+ggtitle("With Daisies")+scale_colour_manual(name="Legend",values=c("brown","black","white", "red"))
plot
plot2 <- ggplot(data.frame(I,Barr,Warr,Uarr,Tarr),aes(I))+geom_line(aes(y=Barr,colour="% Black"))+ geom_line(aes(y=Warr,colour="% White"))+ geom_line(aes(y=Uarr, colour="% Bare Ground"))+geom_line(aes(y=Tarr/25,colour="Temperature"))+ylab("Temperature & Albedo")+xlab("Time")+ggtitle("With Daisies")+scale_colour_manual(name="Legend",values=c("brown","black","white", "red"))
plot2
plot_ly(z=~TEMP)%>% add_surface() %>% layout(title="With Daisies",
scene = list(
xaxis = list(title = "Time"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "T") ))